Edit on GitHub

A Detailed Look ...

... at spatial interpolation, temporal integration, and input/output

  • put together uvetc dictionnary
    • read gridded velocity output (Udata, Vdata)
    • read trajectory output (float_traj*data)
  • interpolate U,V along trajectory from gridded output
    • compare with u,v from float_traj*data
  • compute whole trajectory using OrdinaryDiffEq.jl
    • compare with x(t),y(t) from float_traj*data

Notes: For additional documentation see e.g. 1, 2, 3, 4

1. Import Software

using IndividualDisplacements, OrdinaryDiffEq
p=dirname(pathof(IndividualDisplacements))
include(joinpath(p,"../examples/recipes_plots.jl"))
include(joinpath(p,"../examples/example123.jl"))
include(joinpath(p,"../examples/helper_functions.jl"));

2. Reload Trajectory Output

from MITgcm/pkg/flt

get_flt_ex_if_needed()
dirIn=joinpath(p,"../examples/flt_example/")
prec=Float32
df=read_flt(dirIn,prec) #function exported by IndividualDisplacements
plt=PlotBasic(df,300,100000.0)
0 1×10 5 2×10 5 3×10 5 4×10 5 0 5.0×10 4 1.0×10 5 1.5×10 5 2.0×10 5

3. Read Gridded Variables

via MeshArrays.jl and into a dictionary

uvetc=example2_setup()
Dict{String,Any} with 13 entries:
  "t0"   => 0.0
  "YC"   =>   data type   = Float32…
  "v1"   =>   data type   = Float64…
  "mskW" =>   data type   = Float64…
  "dx"   => 5000.0
  "u1"   =>   data type   = Float64…
  "XG"   =>   data type   = Float32…
  "YG"   =>   data type   = Float32…
  "u0"   =>   data type   = Float64…
  "mskS" =>   data type   = Float64…
  "t1"   => 6.48036e7
  "XC"   =>   data type   = Float32…
  "v0"   =>   data type   = Float64…

4. Visualize Velocity Fields

plt=heatmap(uvetc["mskW"][1,1].*uvetc["u0"][1,1],title="U at the start")

plt=heatmap(uvetc["mskW"][1,1].*uvetc["u1"][1,1]-uvetc["u0"][1,1],title="U end - U start")
10 20 30 40 20 40 60 80 U end - U start - 0.00000003 - 0.00000002 - 0.00000001 0 0.00000001 0.00000002 0.00000003

5. Visualize Trajectories

(select one trajectory)

tmp=df[df.ID .== 200, :]
tmp[1:4,:]

4 rows × 14 columns

IDtimelonlatdepijketaNuVelvVelthetasalttile
Int64Int64Float32Float32Float32Float32Float32Float32Float32Float32Float32Float32Float32Int64
12003600187630.021119.0-1406.2538.0264.723793.00.114312-0.03859630.01014730.29864335.01
220072001.87491e521155.8-1406.2537.99824.731153.00.114318-0.03871940.01028620.29864335.01
3200108001.87351e521193.1-1406.2537.97024.738613.00.114313-0.03884170.0104430.29864335.01
4200144001.87211e521230.9-1406.2537.94224.746193.00.11431-0.03896310.01059940.29864335.01

Super-impose trajectory over velocity field (first for u ...)

x=uvetc["XG"].f[1][:,1]
y=uvetc["YC"].f[1][1,:]
z=transpose(uvetc["mskW"][1].*uvetc["u0"][1,1])
plt=contourf(x,y,z,c=:delta)
plot!(tmp[:,:lon],tmp[:,:lat],c=:red,w=4,leg=false)
0 1×10 5 2×10 5 3×10 5 4×10 5 0 5.0×10 4 1.0×10 5 1.5×10 5 2.0×10 5 - 0.00003 - 0.00002 - 0.00001 0 0.00001 0.00002 0.00003

Super-impose trajectory over velocity field (... then for v)

x=uvetc["XC"].f[1][:,1]
y=uvetc["YG"].f[1][1,:]
z=transpose(uvetc["mskW"][1].*uvetc["v0"][1,1])
plt=contourf(x,y,z,c=:delta)
plot!(tmp[:,:lon],tmp[:,:lat],c=:red,w=4,leg=false)
0 1×10 5 2×10 5 3×10 5 4×10 5 0 5.0×10 4 1.0×10 5 1.5×10 5 2.0×10 5 - 0.00005 - 0.000025 0 0.000025 0.00005

6. Interpolate Velocities From Gridded Fields

uInit=[tmp[1,:lon];tmp[1,:lat]]./uvetc["dx"]
nSteps=Int32(tmp[end,:time]/3600)-2
du=fill(0.0,2);

Visualize and compare with actual grid point values – jumps on the tangential component are expected with linear scheme:

tmpu=fill(0.0,100)
tmpv=fill(0.0,100)
tmpx=fill(0.0,100)
for i=1:100
    tmpx[i]=500.0 *i./uvetc["dx"]
    ⬡(du,[tmpx[i];0.499./uvetc["dx"]],uvetc,0.0)
    tmpu[i]=du[1]
    tmpv[i]=du[2]
end
plt=plot(tmpx,tmpu,label="u (interp)")
plot!(uvetc["XG"].f[1][1:10,1]./uvetc["dx"],uvetc["u0"].f[1][1:10,1],marker=:o,label="u (C-grid)")
plot!(tmpx,tmpv,label="v (interp)")
plot!(uvetc["XG"].f[1][1:10,1]./uvetc["dx"],uvetc["v0"].f[1][1:10,1],marker=:o,label="v (C-grid)")
0.0 2.5 5.0 7.5 10.0 - 4×10 - 5 - 2×10 - 5 0 2×10 - 5 4×10 - 5 u (interp) u (C-grid) v (interp) v (C-grid)

And similarly in the other direction

tmpu=fill(0.0,100)
tmpv=fill(0.0,100)
tmpy=fill(0.0,100)
for i=1:100
    tmpy[i]=500.0 *i./uvetc["dx"]
    ⬡(du,[0.499./uvetc["dx"];tmpy[i]],uvetc,0.0)
    tmpu[i]=du[1]
    tmpv[i]=du[2]
end
plt=plot(tmpx,tmpu,label="u (interp)")
plot!(uvetc["YG"].f[1][1,1:10]./uvetc["dx"],uvetc["u0"].f[1][1,1:10],marker=:o,label="u (C-grid)")
plot!(tmpx,tmpv,label="v (interp)")
plot!(uvetc["YG"].f[1][1,1:10]./uvetc["dx"],uvetc["v0"].f[1][1,1:10],marker=:o,label="v (C-grid)")
0.0 2.5 5.0 7.5 10.0 - 3×10 - 5 - 2×10 - 5 - 1×10 - 5 0 u (interp) u (C-grid) v (interp) v (C-grid)

Compare recomputed velocities with those from pkg/flt

nSteps=2998
tmpu=fill(0.0,nSteps); tmpv=fill(0.0,nSteps);
tmpx=fill(0.0,nSteps); tmpy=fill(0.0,nSteps);
refu=fill(0.0,nSteps); refv=fill(0.0,nSteps);
for i=1:nSteps
    □(du,[tmp[i,:lon],tmp[i,:lat]],tmp,tmp[i,:time])
    refu[i]=du[1]./uvetc["dx"]
    refv[i]=du[2]./uvetc["dx"]
    ⬡(du,[tmp[i,:lon],tmp[i,:lat]]./uvetc["dx"],uvetc,tmp[i,:time])
    tmpu[i]=du[1]
    tmpv[i]=du[2]
end
plt=plot(tmpu,label="u")
plot!(tmpv,label="v")
plot!(refu,label="u (ref)")
plot!(refv,label="v (ref)")
display(plt)
/home/travis/.julia/packages/GR/yMV3y/src/../deps/gr/bin/gksqt: error while loading shared libraries: libQt5Widgets.so.5: cannot open shared object file: No such file or directory
connect: Connection refused
GKS: can't connect to GKS socket application

GKS: Open failed in routine OPEN_WS
GKS: GKS not in proper state. GKS must be either in the state WSOP or WSAC in routine ACTIVATE_WS
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine FILLAREA
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine FILLAREA
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine POLYLINE
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine POLYLINE
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine POLYLINE
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine POLYLINE
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine POLYLINE
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine POLYLINE
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine POLYLINE
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine POLYLINE
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine POLYLINE
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine POLYLINE
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine POLYLINE
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine POLYLINE
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine POLYLINE
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine POLYLINE
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine POLYLINE
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine POLYLINE
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine POLYLINE
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine POLYLINE
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine POLYLINE
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine POLYLINE
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine POLYLINE
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine POLYLINE
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine TEXT
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine TEXT
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine TEXT
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine TEXT
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine TEXT
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine TEXT
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine TEXT
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine TEXT
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine TEXT
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine TEXT
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine TEXT
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine TEXT
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine TEXT
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine TEXT
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine TEXT
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine TEXT
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine TEXT
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine TEXT
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine TEXT
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine TEXT
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine TEXT
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine TEXT
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine TEXT
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine POLYLINE
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine POLYLINE
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine POLYLINE
GKS: GKS not in proper state. GKS must be either in the state WSAC or SGOP in routine POLYLINE

6. Recompute Entire Trajectories

Solve through time using OrdinaryDiffEq.jl with

  • ⬡ is the function computing du/dt
  • uInit is the initial condition u @ tspan[1]
  • tspan is the time interval
  • uvetc are parameters for ⬡
  • Tsit5 is the time-stepping scheme
  • reltol and abstol are tolerance parameters
tspan = (0.0,nSteps*3600.0)
#prob = ODEProblem(□,uInit,tspan,tmp)
prob = ODEProblem(⬡,uInit,tspan,uvetc)
sol = solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
sol[1:4]
retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 4-element Array{Float64,1}:
  0.0               
  0.2007818920948153
  2.2086008130429686
 22.2867900225245   
u: 4-element Array{Array{Float64,1},1}:
 [37.525990625, 4.22379453125]         
 [37.52598925375626, 4.223794946058778]
 [37.52597554129228, 4.223799094139329]
 [37.5258384139952, 4.223840574222589]

Compare recomputed trajectories with those from pkg/flt

ref=transpose([tmp[1:nSteps,:lon] tmp[1:nSteps,:lat]])
maxLon=80*5.e3
maxLat=42*5.e3
show(size(ref))
for i=1:nSteps-1
    ref[1,i+1]-ref[1,i]>maxLon/2 ? ref[1,i+1:end]-=fill(maxLon,(nSteps-i)) : nothing
    ref[1,i+1]-ref[1,i]<-maxLon/2 ? ref[1,i+1:end]+=fill(maxLon,(nSteps-i)) : nothing
    ref[2,i+1]-ref[2,i]>maxLat/2 ? ref[2,i+1:end]-=fill(maxLat,(nSteps-i)) : nothing
    ref[2,i+1]-ref[2,i]<-maxLat/2 ? ref[2,i+1:end]+=fill(maxLat,(nSteps-i)) : nothing
end
ref=ref./uvetc["dx"]

plt=plot(sol[1,:],sol[2,:],linewidth=5,title="Using Recomputed Velocities",
     xaxis="lon",yaxis="lat",label="Julia Solution") # legend=false
plot!(ref[1,:],ref[2,:],lw=3,ls=:dash,label="MITgcm Solution")
25 30 35 40 45 0 10 20 30 Using Recomputed Velocities lon lat Julia Solution MITgcm Solution

This page was generated using Literate.jl.